import netCDF4 
import sys, struct
import math as pm
import numpy as np
import pylab as p
import math
import csv
import matplotlib.pyplot as plt
import h5py
from netcdf_reader import *
from scipy import ndimage
from sklearn.decomposition import PCA
from sklearn.cluster import DBSCAN

INPUT_DATA_DIR = '/home/behollis/DATA/pierre/ocean/'
OUTPUT_DATA_DIR = '/home/behollis/Dropbox/visweek2015/oslines/'

TOTMEMS = 16
X = 152
Y = 152
DIM = 2
 
def writeElapsedTime(start, end):
    import os
    systime = os.times()
    tout = open(OUTPUT_DATA_DIR+'timingDBSCANStirring.{0}.txt'.format(systime[4]), 'w') #sys time
    tout.write('Elapsed time: {0}'.format(end - start))
    tout.close() 

if __name__ == '__main__':
    
    import time
    start = time.time()
    
    output = h5py.File( OUTPUT_DATA_DIR+'stirTermClusBack.hdf5', 'w')
    dclus = output.create_dataset(name='ptclusters', shape=(X,Y), dtype='f')
    
    f = h5py.File(OUTPUT_DATA_DIR +'industrialStirringSlinesTS050.500steps.ss05.2.hdf5', 'r')
    mems = [ d for d in f ]
    
    for x in range( 0, X ): 
        for y in range( 0, Y ):
            tpts = list()
            #tpts = [[],[]]
            print 'calculating terminal clusters for: {0}, {1}'.format(x, y)
    
            for mem in mems:
                dir = mem + '/x' + str(x).zfill(3) + '/y' + str(y).zfill(3)
              
                try:
                    #pt0_0 = f[dir][0][0] #end point in forward integration
                    #pt1_0 = f[dir][1][0]
                    
                    pt0_1 = f[dir][0][-1] #end point in backward integration
                    pt1_1 = f[dir][1][-1]
                    
                    #if math.isnan(pt0_0) is not True and math.isnan(pt1_0) is not True:
                    #    tpts.append([pt0_0, pt1_0])
                    if math.isnan(pt0_1) is not True and math.isnan(pt1_1) is not True:
                        tpts.append([pt0_1, pt1_1])
                    #print [pt0_0, pt1_0]
                    #print [pt0_1, pt1_1]
                except:
                    print 'could not read values from hdf5 file'
                    print '\t for x{0} y{1} {2}'.format( x, y, mem )
                    continue
                
                
            if len(tpts) < 2:
                print 'no pca calculated for land mask,etc. '
                continue  
           
            
            #dcov[x,y] = np.cov(np.array(tpts).T)#pca.get_covariance()
            # let's use eps to be some signficant fraction of
            # the simulation physical domain...we have 152 x 152 cells
            # five percent of the diagonal distance across the domain
            # thus, diag dist = 214.96
            # 0.05 * 214.96 = 10.7
            # NOTE: include this metric and value in the paper
            # min samples is 0.1 of total members, or 2 
            db = DBSCAN(eps=10.7, min_samples = 2).fit( np.asarray(tpts) )
            
            core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
            core_samples_mask[ db.core_sample_indices_] = True
            labels = db.labels_
    
            # Number of clusters in labels, ignoring noise if present.
            n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
            
            dclus[x,y] = n_clusters_
            
            
            
            
    end = time.time()
    writeElapsedTime(start, end)

    f.close()
    output.close()
    
    print 'finished!'